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The Earth's density distribution can be approximately considered piecewise continuous at the 
scale of two-flavor oscillations of typical solar neutrinos, such as the beryllium-7 and boron-8 neu- 
trinos. This quite general assumption appears to be enough to analytically calculate the day-night 
asymmetry factor for such neutrinos. Using the explicit time averaging procedure, we show that, 
within the leading-order approximation, this factor is determined by the electron density within 
about one oscillation length under the detector, namely, in the Earth's crust (and upper mantle 
C*") ■ for high-energy neutrinos). We also evaluate the effect of the inner Earth's structure on the ob- 

served asymmetry and show that it is suppressed and mainly comes from the neutrinos observed 
near the winter and summer solstices. As a result, we arrive at the strict interval constraint on the 
asymmetry, which is valid within quite a wide class of Earth models. 
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I. INTRODUCTION 

The effect of neutrino oscillations in vacuum lies beyond the Standard Model and is thus interesting both from the 
theoretical and experimental points of view. The oscillations in medium have also been studied since Wolfenstein, 
who showed that neutrinos acquire a specific flavor-dependent potential due to coherent forward scattering on the 
tIh ' matter [lj. As a result, the neutrino propagation in medium should demonstrate the conversion from one flavor into 
another (i.e. flavor oscillations), even if the vacuum mixing is negligible. This spectacular result is known as the 
Mikheev-Smirnov- Wolfenstein effect @ and suffices to explain the deficit of observed solar electron neutrinos 
According to Mikheev and Smirnov, the leading-order estimation for the electron neutrino flux depends only on the 
points where the neutrino was born (the core of the Sun) and absorbed (the detector). However, the properties of 
the medium between these two points can also slightly affect the flavor composition of the observed neutrino flux, 
leading, in particular, to the day- night (solar neutrino flavor) asymmetry 0, @. The latter effect being crucial for 
the entire flavor oscillations framework, a number of experiments were set up to catch this slight flavor composition 
variation resulting from the nighttime neutrino propagation through the Earth Some experiments, which should 

be sensitive enough to distinguish the effect, are also planned in the next decade and are now under construction Q. 

However, one should hold in mind that, in order to make a conclusion on the expected presence of the day- 
night asymmetry, one needs not only the experimental data (in terms of event rates and energy spectra) but also a 
convenient theoretical prediction. First experiments in neutrino oscillations were only aimed at outlining the domains 
in the neutrino mass and mixing parameter space which do not contradict the observed rates. For such needs, the 
theoretical estimations obtained using numerical simulations are quite acceptable (see, e.g., [10]). Indeed, although 
numerical experiments can lack perfect accuracy in some regimes and do not give general results, the calculation of 
the desired effect can be performed for the entire parameter space without drastic algorithm changes. Therefore, 
using numerical predictions, one could more or less interpret the experimental data in terms of domains in this 
parameter space which are excluded. But now that various neutrino experiments have yielded a considerable amount 
of data on the vacuum neutrino mixing (Ill - fl3j and its parameters are fixed quite firmly, the question is whether other 
neutrino oscillation effects to be observed are consistent with the framework accepted so far. Namely, one may ask: 
how should one compare the present and future experimental data on the day-night asymmetry with the theory of 
neutrino oscillations in order to make an ultimate conclusion about their agreement within the allowed class of Earth 
and solar models? The latter issue requires a formalism able to give strict constraints (in a form of an interval with 
fixed boundaries) on the predicted asymmetry within certain approximations, but valid within quite a wide class of 
the Earth and solar models. Developing such a formalism constitutes the principal goal of the present paper. In 
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particular, using it, we readily find the constraints on the day-night asymmetry for beryllium and boron neutrinos 
observed in such experiments as SNO, Super-Kamiokande, and Borexino 

Within the neutrino oscillations framework, it is common to use the Schroedinger-like equation to describe the 
spatial variations of the neutrino flavor [l], Within the two-flavor approximation, the Schroedinger problem is 
posed in terms of the 2x2 flavor evolution matrix (operator) R(x,Xq), whose elements are the neutrino flavor 
transition amplitudes after traveling from point xq to x. The evolution equation and the initial condition read, 
respectively, 

dR(x,x ) - , . m 
7^ = -iXA{x)R{x,x ), (1) 

R(x ,x ) = 1, (2) 

where the matrix Hamiltonian is a point-dependent linear combination of the Pauli matrices, 

A(x) = a(x)a 3 + bo\, (3) 
2EV(x) 

a(x) = -cos26» + 6 = sin20 o . (4) 

Am 2 

Here, 9q is the vacuum mixing angle, E is the neutrino energy, Am 2 is the difference between neutrino masses squared, 
and V(x) = ^/2G-pN e {x) is the Wolfenstein potential, which is proportional to the electron density in the medium N e 
and the Fermi constant Gp — 1.17 x 10~ n McV~ 2 . The constant coefficient A = ato 2 /4£J is the reciprocal vacuum 
neutrino oscillation length, up to the factor of 7r, 

AdSC = n = T" • (5) 

Am A 

Equation (p} defines a one-parametric subgroup of SU{2) and hence of 5*0(3), the translation along x being the 
group operation. In this sense, equation (fTJ) is analogous to the Bargmann-Michel-Telegdi equation [ljj in the spinor 
representation 1 1 "i . It is well known that the solution of the matrix linear ordinary differential equation, such as 
([T]), can be represented as a time-ordered exponential (Dyson expansion) [l6j . However, in the case of general N e (x) 
profile, the solution in terms of an analytical function of the matrix argument (i.e., without a symbolic operation such 
as time ordering) appears to be too challenging to find. In the recent investigations, considerable progress was made 
in finding exact solutions of Eq. (fTJ) in certain special cases [13, EH; nevertheless, the general approach to this kind of 
equation still remains to be approximate. 

It was first naturally accepted that numerical simulations could give exhaustive results here and, in a sense, resolve 
the analytical difficulties that arise in the context of such equations. For instance, a numerical approach was employed 
m ll9( to find the domains in the neutrino mixing parameter space where the day-night asymmetry should be 
potentially observable. However, as mentioned above, estimations obtained with numerical techniques are not always 
reliable enough. In particular, in the large- Am 2 regime, which has now proved to be realistic [l2l Il3l|. the oscillation 
length ([5J becomes small, and numerical evaluation of rapidly oscillating solutions introduces large uncertainties. 
The numerical time averaging of the flavor observation rates also becomes inaccurate in this regime. This effect is 
especially strong for low-energy neutrinos, including beryllium neutrinos; probably, due to this fact, the large-Am 2 
and low-energy area was not covered by original numerical simulations [loj . Therefore, obtaining stringent constraints 
on the day-night asymmetry factor for the continuous measurement of small-length neutrino oscillations favors the 
analytical approach. 

Quite a large number of publications are devoted to the analysis of such approximate analytical solutions. Probably 
the most effective technique for finding the approximate solutions of matrix linear differential equations, such as ([IJ, 
is the so-called Magnus expansion [2(J , which is a generalization of the Baker-Campbell-Hausdorff formula [2l[ . This 
approach provides the solution up to any order of approximation, as well as the constraints on the remainder terms 
[22j . Unfortunately, this technique [23T - I261 ] . as well as other general methods (see, e.g., [13, HH), does not readily 
provide a way to find the solution in its explicit and practically usable form, not firmly fixing the Earth model, 
namely, N e {x) density distribution. It is thus desirable to find an approximate analytical solution of equation (J]), 
which is valid under quite general assumptions about the electron density profile N e (x). This idea was developed in 

[US EH. 

As it was mentioned earlier in this section, in our paper, we are not only aiming at finding the relevant approxi- 
mate expressions, but also at analyzing their accuracy and applicability domain. Namely, in section [Hj we find the 
approximate solutions for the flavor evolution matrix inside the Earth, and then, in section Hill we arrive at the ob- 
servation probabilities for the neutrinos of different flavors. These probabilities are finally subjected to the averaging 
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procedure due to the continuous observation of the solar neutrinos throughout the year (Sec. MI B[) . and the results 
of this averaging are discussed in section IIVI The magnitude of the day-night asymmetry appears to be sensitive to 
the non-trivial structure of the Earth's crust under the neutrino detector, so the effect of the crust is paid special 
attention in section IPlI CI The uncertainties of our estimation are also discussed in section HVl in the interesting cases 
of beryllium and boron neutrinos. We analyze the sources of such uncertainties and finally compare our results with 
the numerical simulation. The detailed derivation of the approximate solutions for the evolution operator is moved 
to the Appendix to make the flow of the paper less complicated. 



II. THE DENSITY PROFILE AND THE EVOLUTION MATRIX 



In our investigation, we use the model electron density profile N e (x) with n—1 narrow segments, where it changes 
steeply, separated by n wide though sloping segments. Let us call these segments cliffs and valleys, respectively. Let 
the cliffs be localized near points Xj, j = l,n — 1, namely, occupy segments [x~ ,x~j] of widths ej, where x^ = Xj±ej/2. 

Then the valleys are [Xj_i, xj] and have the widths Lj Xj—Xj—j.. In fact, this kind of model is a good approximation 
for the Earth's densi ty p rofile known in geophysics, where it corresponds to the so-called Preliminary Reference Earth 
Model (PREM) (3^, [H| . In the present section, we will assume that the above segments can be chosen in such a way 
that the inequality 

tj < £ osc < L 3 (6) 

takes place, i.e. the cliffs are narrow and the valleys are wide compared with the oscillation length. In the following, 
we will briefly call such density distribution piecewise continuous. 

Let us note in advance (see Sec. IIVI for details) that for the energies of beryllium neutrinos (E — 0.862MeV, 
£ osc ps 30 km) and, to some extent, boron neutrinos (E ~ 5 — 10 MeV, ^ osc ps 150 — 300 km), the Earth's layers can 
be divided into cliffs and valleys satisfying ^ , as well as into a number of layers whose widths are of the order of the 
oscillation length, while the density variations are small. One can show that the contributions of the latter segments 
can be easily considered in a manner very similar to those of the cliffs, hence we do not consider them until Sec. MI CI 
where they will be paid special attention. In the present section, it is enough to note that the valley-cliff classification 
depends on the neutrino energy. It is also worth saying here that during the night, the neutrino ray traverses different 
paths through the Earth, thus, the lengths Cj and Lj vary. However, assumption (J6j) holds for the most part of the 
night. 

The total flavor evolution operator for our piecewise continuous density profile equals the matrix product of the 
evolution operators for all cliffs and valleys. Within each of these segments, the two small parameters arise: the first 
of them, 



2EV(x) < j 1.0 x IP' 2 , E = 0.862 MeV, 

> 1.2 10 : . E 10 MeV. 



due to the relatively small density of the Earth |29j, |34| , while the second parameter due to the piecewise continuous 
structure of the density profile, 

r _ / 4sc/ 'Lj for jth valley, 
~ I ej/l osc for jth cliff. {S) 



The smallness of parameter 5 will be discussed in Sec. IIVI 

The evolution matrix for each segment, as well as the total evolution matrix, can be subjected to the following 
unitary transformations [25j|: 

R(x,x ) = Z + (x)Y^(x))Ro(x,x a )Y~ 1 (^(xo))Z-(x ), (9) 
Y(ip(x)) = cosip(x) + i<7 3 sinip(x) = exp{i<j3ip(x)}, (11) 



where u)(x) = \J a 2 (x) + b 2 and ip(x) = X j ui(x)dx is the corresponding phase incursion. For the calculations which 
follow, it is also useful to introduce the effective mixing angle in the medium 6{x) £ [0, n/2] [2j, which is defined by 
the expressions 

cos20(x) = -44. sin20(a;) = (12) 

LJ{X) UJ(X) 
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In terms of this angle, 

Z ± (x)=exp{±ia 2 0(x)}. (13) 

The transformation with matrices Z ± (x) locally diagonalizes the Hamiltonian H(x) in the point x. It thus makes 
the complete diagonalization in the homogeneous case N e (x) = const 2]. The transformation with the operator 
Y(ip(x)) isolates the effect of the medium inhomogeneity, the transformed evolution matrix Ro(x,Xo) satisfying the 
equation 

dRo( f> Xo) = -W(x)a 2 e 2i ^R Q (x, x ), (14) 
ox 



where the dot over denotes the gradient 



Due to the fact that the neutrino detector is homogeneous (0(x) = 0), Eq. (fT4|) has a well-defined and physically 
relevant x — > +oo limit for any fixed xq. Moreover, asymptotically convergent behavior of such systems of differential 
equations is stated by the Levinson's theorem [l6[ . 

In the homogeneous case, the equation above is trivial, Rq(x,xq) = 1. However, in the valley [Xj,xJ +1 ], the slow 
change of the density N e (x) enables us to use the so-called adiabatic approximation leading to the same result Q 

R Q (xj +1 ,X+) = 1 + 0(A?7 (5 vancy ), (16) 

where the remainder term is a (generally speaking, non-diagonal) matrix and A77 is the variation of the density 
parameter 77 over the valley. On the other hand, if the Wolfenstein potential undergoes a considerable change within 
the narrow cliff [xj ,Xj] with the phase incursion Aip -C 2ir, then we get 

Ro(x+,xJ) = exp {-ia 2 A0 je 2i ^ x ^} + 0( V S cliS ), (17) 

where A0j = 0(x~j) — 0{xj) = 0(rf) is the jump of the effective mixing angle on the jth cliff. Moreover, one can show 
that within the more accurate 0(r)5) approximation, the above expressions take the form (see the Appendix) 

Ro(xJ +1 ,x+) = exp j-^g [e 2i ^ x ^0(x- +1 ) - e 2 ^ {x ^0{x+)] j + (a V 6 2 ) (valley), (18) 
R (x+,xj) = exp{(-i ( x 2 A^+i(7 1 ^)e 2i,T3,/ ' (x 3 7) }+ 0(rjS 2 ) (cliff), (19) 

where 



/i, =2Xj(y- x-)0{y)Ay = 0( V 5). (20) 

Now let us write the evolution operator for the whole neutrino path. The neutrinos observed during the day are 

created in the point xq inside the solar core, then travel to the Earth, enter the detector in the point x\ and are finally 
absorbed in the point x* inside it. In the nighttime, however, after reaching the Earth in the point x\, the neutrinos 
pass through a number of Earth's layers (valleys) discussed above, and only after that they enter the detector in the 
point x n and are absorbed in x* . Crossing the Sun-to- vacuum interface, as well as traveling inside the Sun, does not 
involve steep electron density changes, moreover, the neutrino propagation is highly adiabatic there (see Sec. IIVI) . 
thus we can treat the whole segment [xo,a;i] as a single valley and use the adiabatic approximation (|16p . As it was 
mentioned earlier, the flavor evolution operator for the whole neutrino path is a matrix product of the evolution 
operators for each segment (each valley and cliff). By the substitution of the approximate solutions (fT5)) and p^|) 
into representation @, after some transformations we find the total evolution operator in the form (compare with 

Eq. mm ) 

R(x*,x ) = ^(z^+y^e^i^-e^ 

X ... X ei^A^g-i^A^giatMxgiasA^tg-iffsesun + 0(^2) (21) 
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Here, the subscript 'Sun' refers to the point xq inside the solar core, where the neutrino is created, and the evolution 
operator inside the neutrino detector is denoted i?dct- The factor n in the remainder term indicates that it contains 
the sum over all cliffs and valleys. Moreover, we use the following notation: 



(xt), = - %T), i = l,n-l, (22) 

(xj), ±9j = d{x+)-d{xj), i=V^T, (23) 



fij = H + -zr= I I 2A (^ - *JM*) + -rrr ) j = i,n-i, (24) 



^j = i J (x j )-ip(x j _ 1 ) = X / uj(x)dx, j = l,n. (25) 



It is also convenient to append definition ((24]) with 

fi n = fi n -0-/2X. (26) 

If the boundary between the Earth's crust and the detector is abrupt, x+ — x~ <C £ sc, then fj, n vanishes. Quite 
analogously, /ii vanishes for the abrupt vacuum-to-Earth boundary. 

By projecting the neutrino state onto the flavor eigenstates, we arrive at the observation probabilities for the 
electron/muon neutrino 

P e , M = | pfe^fy } = l±\sp{R(x*,x )a 3 R\x*,x )a 3 } = (27) 
where, for nighttime neutrinos, 



T = T Dight = i Sp{<r 3 Rl LCt (x*,x+)a 3 R dct {x*,x+ 



x e -io- 2 A6l„„2 e io-i/i„_2 e iCT3Ai/>„„2 e -icr2Afli e icri/ii e icr3Ai/'i e -2icr26(sun 

x g-i<7 3 Ai/>i e io-i£i e -ia- 2 A0i e ~ i<T 3 ^n-i gici/in-i e -io- 2 A0„_i giaijUn g-io- 3 AV>n gif20^ 1 ^8) 

while for daytime neutrinos we have 

T = T day = ^Sv{a 3 R\ ct {x\xi)o 3 R dct {x\xiy^ e °e 1 °^ (29) 

The mixing angle immediately before the detector obviously takes the vacuum value 9q hi this case. 

Due to the homogeneity of the detector substance and its smallness compared with the oscillation length, we easily 
find 

Rdet(x*,x+) = Z+ ct e^ A ^Zr ct = (1 + i^A^e-^ 9 ** + 0(S 2 dct ), (30) 

cr 3^det( x *' a; n) cr 3-Rdot(a;*,^) = 1 - 2kri Af/'dot sin 26> dot + (<5 dct ) , (31) 

where the small parameter 5 de t is the ratio of the detector width -L de t to the oscillation length £ osc . The quadratic 
remainder terms can obviously be neglected. 



III. DAY-NIGHT ASYMMETRY 
A. Finding the probabilities 



In order to evaluate the probabilities obtained above, let us first make the averaging over the phase Atpi, which 
corresponds to the neutrino path between the creation point xq and the Earth. The region of the neutrino creation 
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is extremely large compared with the oscillation length, thus, after this averaging, (cos2aV>i) = (sin2A^i) = with 
a high accuracy. Using this fact together with Eq. (|3"TT) . after averaging (|2"5|) we find 



< - i —'-"^e ") = COS2#Sun 



(Tday) = COS 26»o COS 20 Su 



(32) 
(33) 



The latter expression constitutes the famous result of Mikheev and Smirnov [2]. On the other hand, after averaging 
over Atpi for the nighttime neutrinos, we obtain 



1 



^night -s- -cos26»sunSp{(e 2i(T20 ™ - 2icriA-0detsin26» de t)e 



i(Tl/i„ i(T3Al/>„ — io^ASn-l i(Tl/i„ 



= i(T3 Al/l n - 



x gio-3Ai/>2 e — 2io-2Aei e 2io-i/iig— io-3Ai/>2 g — ia-3Ai/>n-i g-ic2 A9„_i gio-i/j„_i g — ia^Ai/v. gieri^„ 1 

where we have made use of the fact that matrices e 1 " 1 ^ 1 and e ~ I<J2A6 ' 1 commute up to a negligible term of the order 

O(A0!/2l). 

For the calculations which follow, we will use the smallness of the jumps a6*i, . . . , A#„_i = 0(?7) and the parameters 
p,i,...,p, n = 0(rjS). Within the linear approximation in the Earth's density parameter 77, this leads to 



T night (A6i, A0 n _!;/ii, . . . ,n n ) = cos20„ cos20 Sun + Afy ■ ^°'^ h j 



3=1 



an 



Mi • 



night 



A0,p=O j = i 



Ae,p=o 



Partial derivatives with respect to the small parameters are 
<9T n i g ht 



9T n j S ht 



dfij 

In the above expressions, 



A9,p.=0 



A9,p,=0 



-icos20 Sun Sp{(<7 2 e 2iCT3( >" + 2 CT3 AV>detsin20 det )e- 2iff3A ^O 
2 cos 26'sun {sin 26~ cos 2Aip n j — 2A^det sin 20det sin 2Aip n> j } 



cos26>s„„ Sp{(ie iICT2t '" cti + 2A^ d et sin20 dct )e 



-2icr 3 Ai> 71 



'} 



= 2 cos 26'sun {sin 20 n sin 2Aip n j + 2A^det sin 20 det cos 2Atp n j } . 



(35) 



(36) 



(37) 



aVv, 



^>(a; n ) - ip{x, ) = A / cj(x)dx = A-L„,.,(l + 0{rj)), 



(38) 



where L„j = x n — Xj is the distance between the boundary of the jth crossed Earth's shell and the detector, measured 
along the neutrino ray. Finally, by substituting the derivatives and (|3T|) into Eq. ([33]) and using the fact that 
sin 2a-0„.„ = 0, we arrive at the final result 



n-l 



^night = cos 26> Sun jcos 26 n + 2 sin 26 n ^ (a6j cos 2Aip n j + fij sin 2Atp n j) 

3 = 1 
n-l 

™4Ai/>dct sin 26>dct A0,- sin 2Ai() n j j , 



(39) 



which is valid up to the terms of the order 0(r]S 2 ) and quadratic terms 0{j] 2 ). The term including the product of 
the detector width A-0det and the oscillating factor flj cos2ai/v ij is of the order 0(j]S5det) and is thus omitted. 

The above expression provides a generalization of the main result of paper [30j for the case of nonzero-thickness 
cliffs and detector. It should be stressed, however, that Eq. (f3T?)) gives poor information on the effects to be measured. 
Indeed, the neutrino experiments last for years, and thus, Eq. (|39|) may only acquire a predictive power after some 
kind of averaging. The averaging procedure should take into account the axial rotation of the Earth (involving the 
integration over the nights), as well as its orbital motion around the Sun. 
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B. Averaging the probabilities 

The averaging procedure can be performed analytically, if the oscillation phase incursions Aip n j vary by much more 
than 2n during the night. Namely, in this case, one can employ the stationary phase technique (see, e.g., [HI). In 
the case of the beryllium neutrinos (E — 0.862 MeV) traveling through the Earth, the oscillation length is about 
30km, while the depths of the valleys L n j vary by many hundreds of kilometers, so the oscillation phase variations are 
indeed large enough to use the stationary phase approximation. To some extent, the same holds for boron neutrinos. 
However, for both neutrino types, there are layers, to which the stationary phase approximation may not apply. These 
are the Earth's crust immediately under the detector and (for boron neutrinos) the upper mantle. These layers are 
discussed in detail in section ITlI CI and Appendix I A 41 and do not interfere with the picture described in this paragraph. 

Let us consider a neutrino traveling through the Earth, which, according to the PREM model [13], consists of 
a number of concentric spherical shells. The boundary between the valleys Xj corresponds to the point where the 
neutrino crosses one of the interfaces between the Earth's shells; let Tj be the radius of this interface (see Fig. [TJ. 
Further, the distances L n ■ between the detector and the points where the neutrino enters/leaves the interface with 
radius Tj are functions of the 'nadir angle' 6n £ [0,n] defined as the angle between the direction to the Sun and the 
nadir in the point of the detector. In terms of the solar elevation angle 8 S [36], the nadir angle is 0n = B s + tt/2. 
The nadir angle, in turn, is a function of the Earth's axial rotation angle r £ [0, 2tt) ('time of day') and the orbital 
motion angle ^ £ [0, 2tt) ('season'). The dependence of the distances on the nadir angle is easily found to be 

L n,j = L±j(e N ) = r n cos6 N ±\Jr}- ~r\ sin 2 N , (40) 
0N < arcsinrj/r„, (41) 

where the upper/lower signs in (|40l) correspond to the neutrino entering/leaving the interface rj (see Fig. [T]). The 
inequality (l4Tj) ensures that the intersection of the neutrino ray with this interface exists. 

In order to find the night average of the electron/muon neutrino observation probabilities, let us note some properties 
of expressions (|39p and (14TJ1) . First, the number of interfaces crossed by the neutrino is defined via the inequality (|41l) . 
so the number of the valleys and, thus, the number of terms entering the sums in (|39[) are changing during the night. 
Therefore, the night average of Eq. (f39|) contains the sum over all interfaces j, with the average of the jth oscillating 
term involving A-0 nj defined as follows: 

(F(9 N )e 2iA M night = / ne N (r))e 2iA ^ e ^» . (42) 

J AT n ight 
0N (r ) <arcsin rj / r n 

Here, F(9n) is some slowly changing function of the nadir angle and Aright is the total duration of the night in 
terms of the Earth's axial rotation angle r, namely, the length of the segment where 6n(t) < tt/2 (the Sun is below 
the horizon). 

Second, the duration of the night is, in turn, a function of the season ^. On the equinox, e.g., AT n i g ht = tt, while on 
the winter solstice, Ar n i g ht — > max. However, the summer nights are just as long as the opposite winter days, so that 

AT night (<T + 7r) = 2lT - AT n i ght (?) , (43) 

and the total duration of the nights over all the year is exactly half the year. Therefore, the averaging over the year 
of days should involve the division by the total duration of the nights, i.e., 7r7V s . For A s 3> 1, the summation over 
the nights can be replaced by the integration, 

? 

and the averaging formula for the terms containing the phase incursion Aip n j takes the form 

(F(e N )e 2iA ^-) nig ht, y car = ^ / d^ J &T F(6 N (t, ?)) e 2iA ^ . (45) 

6n (t, 1 ?) <arcsin rj /r n 

Now we are able to apply the stationary phase technique to such an integral containing the rapidly oscillating 
exponential. Indeed, let us use the expression, which is valid for smooth functions f(x) and S(x) defined on segment 
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FIG. 1: Radial distribution of the electron density N e (r) inside the Earth and the neutrino path through it. The figure 
demonstrates the cross section of the Earth which contains the nadir DO, the center of the Earth O, and the neutrino ray 



[a, b] containing a single non-degenerate stationary point xq <E (a, b) such that S'(xq) = 0, S"(xq) ^ |3£ 

^ 27F ff r \ p iAS(*o)+i(7r/4)s g nS"(>o) , /Q/)e' ASfa ) 

X\S"(x )\ I[Xo)e + i\S'(y) 



+ 0(A~ 3 / 2 ), A^+c 



(46) 



The two leading terms come from the stationary point and the boundary, respectively. However, in the application 
to the integral (|42p . the boundary term vanishes. Indeed, the boundary of the integration domain corresponds to the 
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neutrino ray being tangent to the interface with radius rj, hence, d T &il> n j oc d T L nt j(Qj^(r)) — >■ oo, and the boundary 
term is absent. On the other hand, the stationary point is obviously achieved at midnight, when the nadir angle 
On — > min (the Sun is in its lowest position) , so the integration over the night yields 



F(e N (r))e M *--( e »W' — = 

7T 

©n (t) <arcsin r j / r n 



1 



7rA|<9 2 L ni 



2iAj/>„,jT i7r /4 



+ 0((\L ntj )-?), (47) 



midnight 



where the two possible signs before m/A correspond to L n j — ■ (see Eq. (00])). The principal point here is that 
the second derivative d^L n j at midnight is suppressed for the inner Earth's shells, 



[9^ nj (0 N (r))] midniglit 
"dX n j(0 N )" 



d(cos On) 



= ± 



midnight 



f&( r\ n din j (On) 

d; cos N ' t 
d(cosON) 

L n ,j{®¥l) 

r V r n- sm2 Or 



midnight 



midnight 



(48) 
(49) 



Now let us use the expression of the nadir ang le On via the Earth's axial tilt e — 23.5°, the latitude of the detector 
X <= [— tt/2, tt/2], and the season ? S [0, 27r] [36| . 



cos On (t, <;) = cos x sin ? sin r + cos e cos x cos ^ cos t + sin e sin x cos 



(50) 



where = corresponds to the winter solstice in the northern hemisphere. The minimum values of On are achieved 
at midnights corresponding to r = r m idnight, 



tanr m idnight(0 



tan^ 
cose ' 



ddnightM = sgn(cos?) 



cose 



V cos 2 e + tan 2 <^ 



Using these expressions, we find the derivative of cos On at midnight 

[d 2 (cosO N )] mi 



cos x cos 2 e + sin 2 <j sin 2 e 



Inudnight | cos ^ | ^ ^ " + ^2 ~ 

Finally, the integral over the night (|47l) takes the form 

dr 



EE -Af(,). 



6n (t) <arcsin rj/r n 



1 



(r 2 /r 2 - sin 2 On) 1 / 4 
V 7rAL nj -(0 N ) 



F(e N )e 



2iAi/i„ i; ,=Fi7r/4 



(51) 



(52) 



(53) 



midnight 



On the other hand, the midnight stationary (minimum) values of the nadir angle ©N^midnight) vary throughout the 
year (see Eq. (|5"TjO . being the smallest on the winter solstice (the darkest midnight) and the largest on the opposite 
summer solstice (the lightest midnight). Therefore, the right side of Eq. (|53"|) is still containing a rapidly oscillating 
function of the season <r, and we can perform another isolation of the stationary points, namely, of the two solstices 
<T = 0,7r: 



(F(O n ) 



night, year 



d£ 
2tt 



^(© N (r,^))e 2iA ^ e ^» — 

7T 



— V - 

2 ^ ^ v / ^(0|a 2 cosO N (r midnight (?))l 
<9 2 cos N (r m idnight (?)) 

s ee sgn{a 2 cosO N (T midn i ght (<;))} = 



0fo 



r„ sin 



in N ) 



@N(T,^)<arcsinrj /r™ 



r 2 / r 2 — sin 2 O r 



AZ/„j(©n) 



-F(0N)e 



2iAi/>„,j+is'(s-l)7r/4 



midnight 



sin(e — x) tane, ^ = (winter solstice), 
sin(e + x) tane, <; = n (summer solstice), 



1, ?-o, 
l, <; = 7r, 

s' = sgn{i nj (0 N ) - r n cos N } = ±1 for L n ,j = L± p 

'l> £>o, 

o, e<o. 



0(0 



(54) 

(55) 

(56) 
(57) 
(58) 
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In the above expressions, the Heaviside theta function $(rj — r n sinOKr) ensures that the stationary point for the jth 
interface is really reached, while the additional sign s' is explicitly introduced to avoid zb-expressions originating from 
Eq. (|40|) . The signs specified in Eq. (|54|) are valid in the northern non-tropical and non-polar latitudes x € ( £ > 7r/2 — e). 
Indeed, the typical neutrino detectors (SNO, Borexino, Super-Kamiokande) are situated in temperate latitudes. Here, 
the midnight solstice nadir angle is ON(7" m idnight) = X + se an d the prefactor j\f(s) = cose cos x- With the use of this 
fact together with the above general averaging formula, one finally arrives at the year average of the night transition 
probability flUJ), 

n-l 

(7 1 night)night,ycar ~ COS 26» Su n COS 26~ + 2 COS 26 S un Sm 26~ ^ A °J 

3 = 1 

0(r,- - r„ sin(x + se)) Jrf /r 2 n - sm 2 { X + se) 
^ ' ' ) \ ■cos{2aV'„ j +s'(s- 1)tt/4}, (59) 



s f^ 1 27r^/sin e cos x sm(x + se) AL nj (x + se) 
s' = sgn{L„j(x + se) -r„cos(x + se)}. (60) 

In the above expression, we have omitted the terms resulting from averaging the (r/S) terms in Eq. (|39[) , since they 
are extremely small. The phase incursions Atp n j should obviously be taken at 0n = X + se, i.e. at solstice midnights. 
For the Borexino detector situated in the Gran Sasso laboratory, with 6 — +42.5°, the prefactor in (I5T)1) which does 
not depend on j amounts to be 

winter solstice (s = —1), 

(61) 




2-7T -^/sin e cos x sin(x + se) 1 0.31, summer solstice (s = +1 



For the Super-Kamiokande detector, 9 = +36.2°, and the prefactor equals 0.60 and 0.30 for the two stationary points, 
respectively. 

Let us briefly note that, in the tropical latitudes |x| < £, additional stationary points appear; in particular, on the 
Equator x = 0, they correspond to the equinoxes. Two additional points meet on the winter solstice, when x - > +c 
(on the Tropic), and one encounters a degenerate stationary point [35] . Although we have found the analytical 
expressions for the year averages in the tropical and equatorial zones (|x| < e), we do not present them here due 
to their mathematical complexity and to the fact that the actual neutrino detectors are situated in the temperate 
latitudes. We confine ourselves to saying that, from (|61[) . one can infer the amplification of the winter solstice 
contribution, as one approaches the Tropic. 

One should also mention that, for low-energy neutrinos (E < 1 MeV), in certain seasons, the phase incursions 
Aip n j for a certain layer j may differ by approximately a multiple of 2n on the successive nights. As a result, the 
contributions of these nights will not cancel each other, quite similarly to those of the nights near the solstices. This 
may be called a parametric resonance and, in principle, will lead to local anomalies in the observed neutrino flux, but, 
as one may see from Sec. IIV1 the day-night asymmetry for low-energy neutrinos is quite small and its anomalies are 
even more challenging to observe. Thus, we do not pay these additional effects much attention here. 

Finally, we are left with the following conclusion. The terms entering Eq. (|39p . which contain the oscillating 

functions of the phase incursions 2Aip n j, are suppressed as O ( jr-jj^ — = O {j J ~'^j after averaging over the year; 

within the leading approximation, the resulting averages (|59|) come from the stationary phase points achieved on the 
winter and the summer solstices. The suppression becomes stronger for the inner Earth's shells. 

It is spectacular that all the terms of the order 0(r]6) in the expression (|3"9")l . including the one corresponding to 
the detector, become 0(r]5 2 rj /r n ) after the averaging. The terms of the order 0(rj), which are proportional to A0j, 
become 0(r]Srj/r n ), respectively. Finally, we are left with the average value 

. r 



(T night ) = cos 26 Sua cos 29- + O Irvfi-L . (62) 

By substituting this result together with the daytime average value (|33ll into expression (f27|) for the neutrino obser- 
vation probabilities, we arrive at the day-night asymmetry factor 

s 2((P e , night ) - (P, day )) = _ W s^ 20 2EV(,-) f rA 

(Pe ,night) + (-Pe.day) 1 + ?MSW COS 29q Am \ T n J 

where Tmsw = cos 2#o cos 2#s U n = a vac asiin/wsun defines the observation probabilities for the solar neutrinos due to 
the Mikheev-Smirnov-Wolfenstein effect [2| and V(x~) is the Wolfenstein potential in the Earth under the detector. 
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One should hold in mind that the asymmetry factor defined as (|63j) may not be directly measurable in the neutrino 
experiments, depending on the detection mechanism. Indeed, definition (|63p differs from the one preferred by the 
experimentalists, the latter being 

^(exp) _ 2(N Dight - TVday) ^ 
-^night ~t~ ^Xday 

where iVday.night is the number of neutrino events observed during the day/night. These numbers may not correspond 
to the electron neutrino fluxes. For example, scattering experiments, which cannot separate the charged and the 
neutral current interactions, are unable to directly measure the electron neutrino flux. To compare prediction (|63p 
with such experiments, one should reinterpret the observed event rates -/Vday,night in terms of the electron neutrino 
fluxes using some theoretical assumptions. 

Let us provide an example of the connection between the day- night asymmetry factors defined as and (|64[) . 
Namely, in the case of the neutrino-electron scattering experiment, such as Borexino, the incident neutrinos produce 
recoil electrons inside the detector, and the scattering cross sections for such processes are well-known [13,|3l], including 
one-loop corrections [39| . The resulting ratio of the total electron/muon neutrino detection cross sections is a function 
of the neutrino energy E (as well as of the threshold T m - m of the recoil electron detection). For monochromatic 
beryllium neutrinos and the actual Borexino's threshold, this ratio is pj| 

a(v e ;E)/a(^;E) ^ 4.5, E = 0.862 MeV, (65) 

thus, the ratio of the event rates is 

Nday _ (-Pe,day)0'(l / e) + (1 - (-Pe,day))0'(^) 



knight (-Pe.nightM^e) + (1 - (•P e ,night)M I/ l*) ' 

The 'experimental' day-night asymmetry factor is then easily found to be 



(66) 



,( exp ) _ . {°j"e) -g(ty))(l +T MS w) ( „ 7] 

dn ~ dn K^ e )-(7(^))(l + rMSw) + 2a(^) [0<) 



and, for beryllium neutrinos (Tmsw ~ 0.09), we obtain 

A dn P) ~ 0.66A dn (E = 0.862 MeV). (68) 

For boron and other types of neutrinos, which have continuous energy spectrum, the expression for the 'experimental' 
day-night asymmetry factor involves the integration over the neutrino energy, 



l( cx P ) = J p(E)dE x A dn (E) {o(y e ; E) - a(^;E))(l + T MS w(E)) 
dn / P {E)dE x {(a(u e ; E) - a{v^E)){\ + T MSW (E)) + 2a{v fl ; E)} ' 



Here, the cross sections <j{v e ^\ E), the 'physical' asymmetry factor A dn , and the Mikheev-Smirnov-Wolfenstein factor 
Tmsw depend on the neutrino energy, and p(E) is the normalized energy distribution of incident neutrinos. As one 
can see, the integration can be easily undertaken numerically using our asymmetry prediction (1631) and the expressions 
for the effective total cross sections, if the recoil electron detection threshold T m i n is known (see, e.g., [38(). 

For those experiments which directly observe the charged-current electron neutrino events, one formally sets 
ct(^ m ; E) — > 0. Then one has A^* v ^ — A dn for monochromatic neutrinos and 



,(ex P ) _ J p(E)a(is e ; E)dE x A dn (E) (1 + T MSW (E)) 
J p(E)a(u e ;E)dE x (1 + T MSW (^)) 



A (cx P) = j ^^V^^^dnv^y^^^M^;; (charged current only) (70) 



for continuous-spectrum neutrinos. 



C. The effect of the crust 



The estimation ([63| shown above is substantially based on the piecewise continuous structure of the density profile 
and shows that the asymmetry should depend only on the density of rock immediately under the detector, i.e. in 
the Earth's crust. At the same time, for beryllium neutrinos, the actual width of the crust is comparable with the 
oscillation length, and neither the valley nor the cliff approximation is valid for this layer. For boron neutrinos, the 
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crust can be considered a cliff, however, the oscillation length becomes comparable with the thickness of the upper 
mantle, as well as with thicknesses of transition zones (see Fig. [T] and [3J|). Moreover, the phase incursions are 
small within the crust (and the upper mantle), and so are their time variations, hence, we are unable to use the 
stationary phase approximation. However, we are still able to account for the effect of these near-surface layers on 
the observed day-night asymmetry factor (|63l) . relying upon relatively small density variation Arj within them. This 
feature, together with the bounded layers' thickness, makes it possible to find the closed form of the approximate 
flavor evolution operator for the crust (the upper mantle) [x^_ 1 , x~] (see Appendix IA 4|) . 



Roix-^t^) = avU-iosfi + tota^^-Jy+Otf), (71) 

X~ 

(3 + ia = { e(y)e 2iWv) -^ x "-^dy, a, (3 = 0(at}) e R. (72) 



This result formally repeats the cliff approximation ([T9")) up to the substitution A0j — > f3, /ij — > a. Using such a 
substitution, the transition zones, whose widths are comparable with the boron neutrino oscillation lengths, can be 
safely replaced by the cliffs, preserving the form of expression (|39p . Finally, to account for the effect of the crust (and 
the upper mantle), as well as the effects of the transition zones, we should make the following modification in (|39[) : 



n-l 



Tnight = cos 26> Sun jcos 26 n + 2 sin 26 n ^ (a8j cos 2ai/>„ j + fij sin 2Aip n j) 

3 = 1 
n-l 

-4Ai/>dct sin 26> dct A0j sin 2ai/>« j } + AT night , (73) 

3=1 

AT n i g ht = cos26' Sun {2sin26' ri ;(/?cos2A?/;„ in _i +asin2Ai/) ni „_i) 

-4A^dct sin26'dct(/3sin2A'(/) ni „_i - a cos 2a^/v„-i)} 

X n X n 

= 2 cos 26» Sun sin 26~ J 6(y) cos 2Aip(y) dy - 4 cos 26» Sun A'0dct sin 26» dc t / 6(y) sin 2 Aip(y) dy, (74) 

where Aip(y) = ip(x~) — ip{y)- The leading 0(rj) correction of the crust (the upper mantle) to the day- night asymmetry 
factor (ffi3")) reads 

. ?msw sin 2 26» 2E f ■ 

AA dn = —— — — — — / V (y cos2a^ (y )dy. 75 

1 + J msw cos zt/ &m J 

Again, this result should be subjected to the averaging procedure to acquire the predictive power. However, due to 
the boundedness of the cosine, the correction is easily estimated (both before and after averaging), 



(A^dn) 



^4d„ 



/ |*(y)|d» ~ (76) 

V{x n ) J V(x n ) 



„+ 



IV. DISCUSSION 



Let us begin with the review of the approximations we were using in the above calculations of the day-night 
asymmetry. First, the Sun was considered one big valley, i.e. the adiabatic approximation was used for it. The 
non-adiabatic corrections are of the order of the adiabaticity parameter 7 = \0(x)\/Xu> @. Using the fact that the 
typical spatial scale of the solar density variation is associated with the solar core radius Ro ~ 0.1i?s un ~ 7 x 10 km 
j40(, we find 

<JlO- 5 , £J~lMeV, 
7Sun ~ I 5 x 10- 4 , £;^10MeV. 1 ' 
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Thus, these corrections can be safely neglected. 

Neglecting the finite width of the detector introduces a relative error of the order O ^-j^ 5 J , compared with the 

leading term (|63l) . as one can see from (|39l) . This correction is minuscule even for beryllium neutrinos (^ sc ~ 30 km), 
since the detector sizes are now Ldet < 1 km. Moreover, this type of correction is subjected to additional suppression 
due to time averaging. 

The valley-crust approximation considered in Sec. Hi] relies upon the smallness of parameters rj and S. The relative 
error of the linear approximation in rj is of the order of 0(rf), hence, this approximation works fine for beryllium 
neutrinos and quite well for boron neutrinos (see Eq. ([7]); note that the maximum values specified for rj correspond 
to the inner core, while for typical neutrino detector latitudes, the Sun never descends low enough to shine through 
it). Thus, for boron neutrinos, the error of the linear approximation is within 10%. Further, the error due to the 
finite width of the jth valley is estimated as O(ArjjSj) — 0(Arjj £ osc /Lj), where AT)j is the density variation on the 
jth valley and Lj is the width of the valley. On the other hand, if the width of the kth cliff ^osc, the resulting 

error is 0(rj8k) = 0(rj Lk/l osc ). Within the PREM model, the cliffs are abrupt (Sf. = 0), so we will not discuss the 
numerical values of the corresponding errors. 

If there are deep layers in the PREM density distribution, whose widths are of the order of the oscillation length, 
they can be considered neither valleys nor cliffs. However, if the density variations Arj over these layers are small, we 
can use the approximation considered in Sec. IIII Cl and Appendix IA 41 which formally reduces these layers to the cliffs 
(see Eq. (fTTj) ). The error of such an approach is measured by the magnitudes of parameters a, (3 in (fTTj) . i.e. the error 
is of the order of Arj. The contributions of such layers are additionally suppressed after averaging, due to their depth 
(see below). In the case of the Earth, there is a layer of this type. Namely, the core-mantle transition zone is about 
200 km wide [Hj], which is of the order of the typical boron neutrino oscillation length. However, the corresponding 
density variation lies within only A 77/77 ~ 5 — 7%, and the time-averaged contribution is obviously negligible. For 
beryllium neutrinos, the core-mantle transition zone can be considered a valley. 

The contributions from the Earth's layers, which lie much more than the oscillation length under the detector, 
are suppressed due to time averaging, the principal nonvanishing contributions being provided by the winter/summer 

solstice stationary points. These contributions have the relative magnitude O C^-^-j^^j compared with the leading 

term ([63]) (see Eqs. (|39|) and (l59l) ). This suppression is quite strong for dense though deep inner Earth's shells, 
including its core. 

Finally, for the layers which lie within O(£ osc ) under the detector and have the widths comparable with the 
oscillation length, the averaging procedure cannot be performed so as to result in the expression (|59p. Here, we are 
only able to follow the approach of Sec. (jlll CI) , which results in an unaveraged correction (1751) . The latter correction, 
in principle, can be explicitly averaged using numerical integration, which in this case involves only the near-surface 
Earth's structure. The strict constraint on the uncertainty introduced by neglecting the effect of this structure is given 
in Eq. (f?6")) . This expression predicts the 20% uncertainty for beryllium and ~ 30% for boron neutrinos. However, it 
is quite obvious that the time average of the cosine in ([75| is sensitive to the Wolfcnstein potential in the layers within 
O(£ osc ) under the detector; moreover, ([5TJf shows that this sensitivity asymptotically falls off as the inverse depth of 
the layer L n j. Thus, if the resulting expression (p3"|) is used for boron neutrinos, one should substitute for V(x~) 
the Wolfenstein potential averaged over the crust and the uppermost layers of the mantle, using some decreasing 
weighting function. The resulting potential will be 15 — 20% larger than that immediately under the Earth's surface. 
For beryllium neutrinos, it suffices to substitute the Wolfenstein potential in the crust. 

Thus, we may conclude that, according to (j6"3")l and (|75[). the asymmetry has the order 0(rj) and is determined 
by the rock density in the layer of width about £ osc under the detector, i.e. in the Earth's crust (as well as in the 
uppermost part of the mantle for boron neutrinos) . The principal uncertainty of expression (1631) comes from the fact 
that the stationary point approximation we used for the analytical time averaging of Eq. (131))) is inapplicable to the 
layers which lie within several neutrino oscillation lengths under the detector. Another correction comes from the 
time averaging of oscillating contributions of the deep layers; using (I59|) . we can estimate its magnitude as < 3% for 
beryllium neutrinos and < 10 — 15% for boron neutrinos with E — 10 MeV. All other corrections, taken together, do 
not exceed 10%. 

Using the recent data from SNO, KamLAND, and Borexino collaborations [lTl-[T3j|. namely, tan 2 9q sa 0.46 and 
Am 2 s» 7.6 x 10~ 5 eV 2 , and the typical electron densities in the Earth's crust N e / CIUSt ^ = 1.3mol/cm 3 [33[ and in the 
solar core -ZV e (g un ) ~ 100 mol/cm 3 (40| . we arrive at the numerical estimation for the day-night asymmetry factor for 
solar beryllium-7 neutrinos (E = 0.862 MeV) 

A dn ( 7 Bc) = (-4.0 ± 0.9) x 10~ 4 , A^ p) ( 7 Be) = (-2.6 ± 0.6) x 10~ 4 . (78) 

The uncertainty corresponds to the effect of the Earth's crust, which, as mentioned, can be explicitly evaluated by 
numerical averaging of Eq. (|75"|) . For boron-8 neutrinos with E — 10 MeV, substitution of the crust density into (|6"3")l 
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leads to the asymmetry estimation 

^4dn( 8 B) « (2.9 ±0.8) x 10~ 2 (£ = 10 MeV). (79) 

The 'experimental' asymmetry factor for the electron scattering experiment with the recoil kinetic ener gy t hreshold 
Tmin = 4.5 MeV (such as Super-Kamiokandc-III [41]), averaged over boron-8 solar neutrino spectrum [42(, is then 
given by Eq. (j6"9")l . 

^° n xp) ( 8 B) « (1.6 ± 0.5) x 10~ 2 (averaged over energy), (80) 

the fully- numerical calculation yielding A d °* p \ s ~B) = 1.9 x 10~ 2 . We do not present here the asymmetry predictions 
for 13 N, 15 0, and pep neutrinos, since their typical energies are around 1 MeV, while the fluxes at least one order 
smaller than that for beryllium neutrinos jicj . 

If, for boron neutrinos, one substitutes into (|63[) the mean density in the near-surface layers, instead of the crust 
density, the above predictions (|79)) . (|80)l will be about 10% larger (depending on the weighting function chosen for 
the mean density evaluation) but never larger than those obtained by substituting the density of the upper mantle. 
Indeed, Fig. [2)3 shows that the two analytical curves corresponding to the two 'surface densities' discussed here enclose 
the fits obtained using two techniques of numerical simulation. One of these simulations involves the numerical time 
averaging of the leading 0(ri) terms in expression (|39[) (the valley-cliff approximation), while the other one includes 
both the numerical solution of the evolution equation (|14[) and the subsequent time averaging. 

The numerical curves shown in Fig. [2] were calculated for the Gran Sasso laboratory, where the Borexino detector 
is operating; the curves for Kamioka (Super-Kamiokande) and Sudbury (SNO) lie very close to that for Gran Sasso, 
so we do not show them in this figure. Instead, a comparison of numerical results for selected detector latitudes is 
presented in Fig. [3] (again, we do not include the SNO latitude, since the results for it lie very close to those for 
Gran Sasso). We have also included in Fig. [3] the numerical results for the Northern tropic, since, according to our 
analytical estimations, near it, the subleading contributions to the asymmetry may become considerable, which come 
from the solstices (see Eq. ([59])). We will address this issue further in this section. 

From Fig. [3l one may infer that the leading-order analytical estimation ()63|) is in good agreement with the numerical 
results even for boron neutrinos. Contrary to these numerical results, however, our analytical estimation is model- 
independent, in particular, it does not contain the latitude of the neutrino detector. The possible errors of our 
predictions can also be easily estimated. In the case of beryllium neutrinos, the estimation of the errors becomes 
strict enough to result in a fixed-boundary interval (not a confidence interval) for the day-night asymmetry shown 
in Eq. (|78| . which is useful for experimental purposes. Our results are also in agreement with numerical day-night 
asymmetry predictions presented by other authors (see, e.g., fiol]). 

As predicted, the agreement of the numerical results with the analytical one becomes better for low-energy neutrinos. 
We also observe that the asymmetry vanishes for low-energy neutrinos, which is a result of the applicability of 
the adiabatic approximation for such neutrinos. Indeed, within this approximation, the neutrino flavor observation 
probabilities depend only on the creation and absorption points. 

It is worth emphasizing here that such a directly measurable quantity, as the day-night asymmetry factor, does 
depend both on the neutrino regeneration effect inside the Earth and on the Mikheev-Smirnov-Wolfenstein effect 
inside the solar core, where the neutrino is created. Another useful physical quantity, namely, the regeneration factor 
/ rcg , describes solely the Earth effect in the day-night asymmetry [30j . 

. 2 cos 20 Sun 2 cos 26> S un f rsi s 

^dn — . (rr .rp \/ /rog~ -. rp /rcg, (ol) 

1 + (J night + 1-day)/* 1 + 1 MSW 

/ reg « -sin 29 w . (82) 



Indeed, one observes that, unlike the regeneration factor (1821 . the day-night asymmetry (|63[) depends on the solar 
effect manifested in the energy-dependent quantity Tmsw- As a result, at the energy E ~ 2.0 MeV, which corresponds 
to the Mikheev-Smirnov resonance in the solar core, we have cos2#sun = and Tmsw — 0, and, as a consequence, 
the day- night asymmetry factor fl6"3"]) changes sign (see Fig. [5^); on the other hand, the regeneration factor is always 
positive. Thus, it is the resonance inside the Sun which makes the asymmetry, observed on the Earth, vanish. 

From Fig. [5k., [3^, one may observe that beryllium neutrinos (E — 0.862 MeV) would be indeed quite useful for 
the study of the matter effects in the neutrino oscillations, since they correspond to almost maximum asymmetry 
magnitude in the unusual domain of its negativity. Moreover, solar beryllium neutrinos are highly monochromatic (in 
contrast, e.g., to the boron neutrinos) and their flux is considerably larger [4(j. However, we are able to conclude that 
the day-night effect needs at least a 10-20 times improvement of current experimental resolution to be distinguished at 
a considerable confidence level for such neutrinos. In particular, the result of the Borexino experiment ( 7 Be) — 
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(l ± 12(stat.) ± 7(syst.)) x 10 -3 was reported in April, 2011 @ and demonstrates the strong dominance of the 
uncertainties over the expected effect. One may still hope that the 50 kton LENA detector which should come 
into operation around 2020 [|| and is expected to observe as much as 10 4 solar beryllium neutrino events per day, 
could distinguish the day-night effect for such neutrinos. Strictly speaking, the Poisson statistics results in the 
relative errors AN/N ~ 1/y/W, and, for a year- long experiment, one reaches the statistical error of the order of 
1/V365 • 10 4 - 0.5 x 10-3. However, one could employ the adaptive processing of the experimental data, taking into 
account the expected form of the curve T n ight(0N) (see (|39)) and (|40|) ). i.e. the dependence of the asymmetry on 
the nadir angle. This processing technique may be efficient in extracting the day-night effect from under the noise 
even for small event rates. Quite a similar technique was recently suggested in |43[ as a search tool for periodic time 
variations in the 7 Be solar neutrino flux observed at LENA. 

In Fig.[3jD, one can also observe the role of the stationary phase points in the time average of the day- night asymmetry 
factor. Namely, the numerical curve corresponding to the Northern tropic demonstrates specific oscillations which 
come from the amplified contribution of the winter solstice stationary point to the year average of the asymmetry 
factor (see Eq. (f5T?)) ). The curve for the Kamioka mine, which is situated about 1.5 times closer to the Tropic than 
Gran Sasso, also demonstrates oscillations, compared to the Gran Sasso curve. In view of this effect, it would be quite 
prospective to build a detector close the Tropic, which could be able to observe high-energy solar neutrinos with an 
energy resolution about 0.5 MeV. A favorable place for such a high-technology project could be, e.g., near Sao Paulo, 
Brazil (latitude \ — — 23°33', i.e. exactly on the Southern tropic!), especially under the potential support of the local 
university. 

One should hold in mind here that the positions of the interference peaks on Fig. [3] substantially depend on the 
radii of the Earth's shells. Nevertheless, the approximate smoothness of the three numerical curves in Fig. [3] indicates 
the (approximate) stability of the day-night asymmetry factor with respect to slight variations of the parameters of 
the PREM model, namely, the radii of the Earth's shells and the density jumps. As seen from Fig. 02 this stability 
becomes stronger for low-energy neutrinos, as well as for the detectors operating far from the Tropic. 

It is also interesting to study the effect of the local Earth's crust inhomogeneities under the detector on the 
observed day-night asymmetry. Such inhomogeneities could be associated, for instance, with oil-bearing horizons. Let 



the inhomogeneity be described by the variation 8Ne (x) of the electron density over the smooth profile N e (x) 

N e (x)=N e (x) + 5N® (x), SNP (x) = for + Sx\], 



(83) 



where the inhomogeneity size 5x\ <C £ osc . Then the contribution of this inhomogeneity to the asymmetry factor is 
given by Eq. (|75|) . 
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where L^> the depth of the inhomogeneity under the detector. One can see that, although the regeneration 

effect is stronger for higher-energy neutrinos, this effect is insensitive to the near-surface local inhomogeneities for 
such neutrinos due to the smallness of the sine and ratio 2TrSxi/£ osc in the above expression. Therefore, exploration of 
the Earth's crust based on the neutrino oscillations would require an extreme improvement of current measurement 
techniques (44| . and its prospects seem obscure in the nearest future. 



V. CONCLUSION 



Let us make a brief summary of our initial goals concerning the analytical approach to the day-night asymmetry 
and the results of our investigation. We have attempted to develop a framework able to give interval constraints on 
the day-night flavor asymmetry, which are independent of the details of the density distribution inside the Earth. Of 
course, some approximations should be made to analytically obtain such general results; in our case, the principal 
assumptions were the relatively small density of the Earth (rj <C 1) and its spherically-symmetric layered structure 
(manifested in the small parameter 6). Although the actual approximation parameters may be not quite small, using 
the framework developed, we can readily estimate the corresponding inaccuracies in our predictions (as done, e.g., in 
Sec.HVl). 
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Our analysis shows that the day-night asymmetry is insensitive to the structure of the deep Earth's layers, including 
its core; we found that this sensitivity falls off as the inverse layer's depth 1/L n j (see (JHU))- The day- night effect 
averaged over time (i.e. the effect directly measured in neutrino experiments) is principally determined by the mean 
electron density of the Earth within 1 — 2 oscillation lengths under the neutrino detector, i.e. in the Earth's crust for 
beryllium neutrinos, as well as in the upper mantle for boron neutrinos. The corresponding leading-order analytical 
curves plotted for these two electron densities are shown in Fig. [2] and [3] (the solid and the dashed line, respectively), 
together with the results of the numerical simulations. One may see that the leading approximation works quite well 
within the energy range E ~ 0.5 — 12 MeV under major interest in the field. Moreover, in Fig. [^p, one may notice 
that high-energy neutrinos, as it was expected, 'sense' the deeper Earth's layers. Indeed, it is indicated by the fact 
that the high-energy segment of the numerical curve approaches the dashed theoretical curve plotted for the (higher) 
upper mantle density. If, in view of this fact, one substitutes into the leading-order expression (|63l) the mean Earth's 
density within 1 — 2 oscillation lengths under the detector, the resulting estimation will agree with the numerical one 
within 10%. Such an accuracy of the day-night effect measurements is yet to be achieved in the future experiments, 
such as LENA 0. 

Further, the theoretical analysis predicts next-to-leading-order corrections to the day-night effect, which may arise 
due to the stationary phase points during the year, at which the nighttime neutrino oscillation phase freezes. These 
points occur on the winter and the summer solstices and are spectacular for the fact that their contributions do 
not vanish after (arbitrarily) long observations. These contributions are very small for low-energy neutrinos, as well 
as in the temperate latitudes, however, in the tropical latitudes, they are quite distinguishable. One can see these 
stationary point contributions in Fig. [5Jd (high-energy segment of Fig. [3]), especially for the Tropical curve (latitude 
X = +23.5°) demonstrating oscillations relative to the other curves plotted for the more temperate latitudes. With 
the present energy resolutions reaching 0.5 MeV [4l|, only the event rates are yet too small to observe this effect. 

It is also worth mentioning that the good agreement of our analytical predictions with the numerical simulations 
in a surprisingly wide range of neutrino energies (i.e. oscillation lengths) is also a byproduct of a number of specific 
features of the actual Earth's density profile. For instance, the large density jumps are lying very deep inside the 
Earth; the crust is quite thin (and, as mentioned, is not a cliff for beryllium neutrinos), however, its density is quite 
low; the layers' widths are incomparable, etc. The framework developed in the present paper, however, is able to 
reveal the situations in which the 'universality' of the prediction (|63[) will not hold; one needs only the general features 
of the density distribution to make such a conclusion. 
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Appendix A: Approximate solutions for the evolution operator 

In this section, we derive the approximate solutions for the flavor evolution operator Rq{x,xq) (see Eqs. @ and 
(|14p) for valleys and cliffs. Our approximations will only deal with the neutrino propagation inside the Earth, since the 
neutrino propagation inside the Sun is highly adiabatic (see the adiabaticity estimations in Sec. IIV|) . i.e., Rq = 1 with 
a great accuracy. Moreover, the Earth regeneration effect under investigation depends only on the Earth's density 
distribution N e {x) and not on the details of the neutrino propagation inside the Sun. 

Now, due to the relatively small Earth's density, which manifests itself as a small parameter rj < 10 _1 , we resort to 
the linear approximation in 77, namely, take the two leading terms of the Dyson series 

( X \ X 

Rq(x,x q ) = Texp J -i^ / 0(y)e 2i ° 3 ^dy\ = 1 - kr 2 / 9(y)e 2i ^^dy + 0(if). (Al) 
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1. Valleys 



In the valleys, we have a smooth and bounded function 8 (y) and a rapidly oscillating matrix exponential e 21173 ^^ = 
cos2ip(y) + i<7 3 sia2tp(y). Then, applying the double integration by parts, we arrive at 
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(A2) 



(A3) 



Let Arj be the total variation of the density parameter rj on the valley. Then the variation of the effective mixing 
angle 6(x) has the order 0(at?), and all the gradients in the above expressions are suppressed by powers of the small 
parameter 6 = £ osc /L = tt/XL, where L is the width of the valley and £ osc is the oscillation length. Namely, 



6 = 0{a9/L) = 0(Ar)/L), 

6 = 0{a6/L 2 ) = 0(a v /L 2 ), 

uj = \Ja 2 + b 2 = 1 + 0(t}), 

ci = ^L§= o(a v /l). 



Using these estimations, one can readily show that 

1 / 



6 CoO 
2\a 3 \ \ bj ui 2 



\\L 2 







[atjI 
XL 2 



0(a v S/L), 



■D 2 J(y) = 0(A V S 2 /L), 
then the matrix norm of the remainder integral in Eq. (|A2 



< I \\V 2 J(y)\\dy = 0(A V 6 2 ). 



(A4) 
(A5) 
(A6) 

(A7) 



(A8) 
(A9) 

(A10) 



On the other hand, the first term in the right side of Eq. (|A2|) consists of two parts, the one proportional to T> y 9(y), 
which is of the order O(Ar) 8 /XL) = 0(atj 8 2 ), and the other proportional to 6(y), which is 0(at] 8). Then, 
summarizing the estimations made, we conclude that 



%) 



a 1\a 3 if){y) 



2ia 3 Xui(y) 



0(Ar] S 2 ). 



(All) 



Moreover, here, due to (jA6|) . we can safely replace ui(y) by unity, and then the approximate solution of the evolution 
equation in the valley takes the form 



Ro(x,x ) = 1 - ^ (e{x)e 2i(73 ^ x) - 9{x )e 2{a3 ^ X0 ^ + 0(a v S 2 ). 
Finally, by neglecting the terms of the order 0((at]) 2 8 2 ), we can write 

R (x,x ) = ex pj-^ (0{x)e 2ia3 ^ -0(x o )e 2iff3V,(xo) )| + 0(at7 S 2 ). 
This expression can also be derived using mathematically strict stationary phase technique [3£ 



(A12) 



(A13) 



2. Cliffs 



On the cliffs, the change of the effective mixing angle is of the order 0(rj), while the spatial scale of this change is 
quite small. Therefore, the change of the phase of oscillations Aip <C 2-7T, and we can expand the exponential in the 
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right side of Eq. (|A1|) in the local phase incursion 

v 

rP{y) - i>{x ) = X f lo(z)Az = X(y - ar ) + 0(r,5), (A14) 



e 2i^(y) _ e 2ia 3 ^(x )( 1 + 2ia 3 A(y - x )) + 0(S 2 ) + 0(r]8). (A15) 
Now, using the boundedness of the total variation of the mixing angle on the cliff, we obtain 

X 

e( y )e 2i ^^dy = e 2ias ^ x ^ [ (1 + 2ia 3 \(y - x a ))9(y)dy + 0(t]S 2 ). (A16) 



This finally leads to the cliff approximation for the evolution operator 

Ro(x,x ) = 1 + (-ia 2 A0 + iCTi/^e 2 ^*") + 0( V 6 2 ) 

= cxp{(-ia 2 A8 + ia 1 ^e 2ia ^ ( - x «')} + 0(? 1 6 2 ) + 0(r] 2 ), (A17) 

where a9 = 9{x) — 9(xq) = 0{rf) and 

X 

» = 2\f(y- x o )0(y)dy = 0(r,S). (A18) 
3. Valley + cliff 

It is also useful to consider a valley [x~j , x~ +1 ] following a cliff [xj , ad"]. Using expressions (| A13I) and (|A17|) for the 
evolution operators on these segments, within the linear accuracy in r], we obtain 

Ro (xj + i,xj) = Rq (xJ +1 ,xf)Ro(xf, x~ ) 



exp 



~ (9(xJ +1 )e 2i<T ^ x i+J - 9(x+)e 2i ^ x t^ + (-i^Aflj + ia lf ij)e 2la ^ > I + 0( V S 2 ) 



= expl-^-6(xj +1 )e 2ia ^ x ^ + (-ia 2 A6j + iaxjij) <F™K*j)\ + O^S 2 ), (A19) 

where A9j = 9(xf) — 9(x~), pj = [ij + 0(x~j)/2\, and \ij is defined analogously to (|A18|) . Now, by substituting the 
above result into representation Q for the evolution operation, we arrive at 



= e i^e(x7 + J e -ia 1 9(x7 + J/2A e ia 3 A^ +1 gia 1 /i 3e -i CT2 Ae 3 . e -i<r 2 0(a:7) + + ( A2 q) 

where Aipj + i = —ip(xj). This expression, in turn, can be generalized to a sequence of n valleys (see Eq. (|2ip), 

4. The Earth's crust 

As one could see from the main flow of the paper, for neutrinos with energies E ~ 1 — 10 MeV, the Earth's crust 
(and, for E ~ 10 MeV, the upper mantle) cannot be considered either a valley or a cliff, since its width is comparable 
with the oscillation length, and the parameter S is of the order of unity. Here, however, another approximation is 
useful, which takes into account the small density variation aO over these layers, as well as their bounded thickness. 
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In terms of parameters r\ and 5, we have 



2EV J 2 x 1CT 3 (beryllium neutrinos, E = 0.862 MeV, the crust), 
Am 2 ~ I 3 x 10~ 2 (boron neutrinos, E = 10 MeV, crust + upper mantle), 



(A21) 



— < 0.3 (both cases), (A22) 
V 

S = — < 5 (both cases). (A23) 

*-osc 

Then, using expansion fjAip , we find the immediate expression for the evolution operator for the crust: 

Ro(x, xo) = t + (ia ia - ia 2 P)e 2i ° 3 ^ + 0( V 2 ) = exp{(-i<7 2 /3 + ia ia )e 2ia ^ x ^} + 0( V 2 ), (A24) 
where real numbers a, /3 = 0(at]) are defined by the expression 

X 

P + ia= [ e(y)e 2i W y) -^ X0 »dy. (A25) 



Note that the form of this approximation (|A24[) coincides with the cliff approximation (|A17|) . up to the coefficient 
substitution A0 j3, fi a. In particular, the cliff approximation is restored in the S — > limit. 
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FIG. 2: Comparison of analytical expression (163[) for the day-night asymmetry factor with the results of numerical simulations 
for Gran Sasso (x = +42.5°). Solid/dashed curves: Analytical estimations (|63|l with the crust/upper mantle electron densities 
substituted; Circles (0) : The result based on numerical averaging of analytical expression (|39p ( 'valley- cliff ' approximation); 
Crosses (x): Fully numerical result (numerical solution of (|14|l and subsequent numerical averaging). Subfigures (a) and (b) 
demonstrate the low- and high-energy segments of the curves, respectively. 
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FIG. 3: Comparison of analytical expression (|63[) for the day-night asymmetry factor with the results of numerical simulations 
for different detector latitudes. Solid/dashed curves: Analytical estimations (|63p with the crust/upper mantle electron densities 
substituted; Crosses (x): numerical result for Gran Sasso (x = +42.5°), Pluses (+): for Kamioka (\ = +36.4°), Bullets (•): 
for the Northern tropic (x = +23.5°). Subfigures (a) and (b) demonstrate the low- and high-energy segments of the curves, 
respectively. 



